Una distribución de probabilidad de una variable aleatoria es una función que asigna a cada suceso la probabilidad de que este suceso ocurra. Se puede comprender, de forma sencilla, como la frecuencia teoríca, o cómo se espera que varíen los resultados.
La distribución de probabilidad más sencilla es la distribución uniforme, que asigna la misma probabilidad a cada suceso de ocurrir.
La distribución uniforme discreta es una distribución uniforme, los sucesos son equiprobables, donde el espacio de probabilidad es finito, es decir, hay un número finito de casos o resultados diferentes.
Si tenemos \(k\) resultados distintos, asignaremos a cada uno la probabilidad \(1/k\), ya que la probabilidad de todos los posibles sucesos tiene que ser 1:
\[\sum_{i=1}^k P[X = x_i] = 1 \]
A los posibles valores que puede tomar una
Algunos ejemplos cotidianos, comúmente usados para estudiar la probabilidad, son el de la moneda, \(k=2\), o el dado, \(k=6\).
Podemos ver de forma empírica, simulando, qué ocurre, con 30 lanzamientos de moneda:
set.seed(42)
muestreo_moneda_peque <- sample(1:2, 30, replace = TRUE)
barplot(table(muestreo_moneda_peque))Y con 10.000
set.seed(42)
muestreo_moneda_grande <- sample(1:2, 10000, replace = TRUE)
barplot(table(muestreo_moneda_grande))O 10.000.000
set.seed(42)
n <- 10000000
muestreo_moneda_grande2 <- sample(1:2, n, replace = TRUE)
barplot(table(muestreo_moneda_grande2))Y para obtener probabilidad, dividimos por el tamaño muestral
La probabilidad de sacar cara (suponiendo que la identificamos con un 1) podemos calcularla así
## [1] 0.4999179
NOTA: La función set.seed() determina una semilla para el generador de números aleatorios, que permite replicar una simulación o muestreo.
EJERCICIO 1
Escoger una semilla aleatoria y estimar, mediante un muestreo con diferentes tamaños muestrales, las probabilidades para un dado de seis caras. Ver cómo varía la estimación. Representarlo mediante un gráfico de barras y de sectores.
Después, para la muestra grande, intenta calcular la probabilidad de sacar entre un 2 y un 5.
La noción de distribución uniforme se puede trasladar, de forma bastante intuitiva, al caso continuo. Esta distribución está definida por los extremos del intervalos en que toma valores, \(a\) y \(b\).
En este caso, sería la integral de la función de densidad en este intervalo la que tendría valor 1:
\[\int_a^b f(x)=1, \quad f(x)=\frac{1}{b-a}\]
El espacio muestral en este caso sería
\[S = \{\text{Todos los números entre a y b}\}\]
Podemos obtener valores de una uniforme (0,1) de la siguiente manera:
Y si generamos lo suficiente, tendremos una buena aproximación de la función de densidad:
n <- 1000000
muestra_unif_grande <- runif(n)
hist(muestra_unif_grande, freq = FALSE)
lines(density(muestra_unif_grande), col = "blue")Es fácil ver que el área es 1, ya que en el gráfico se puede apreciar un cuadrado con ambos lados de longitud 1.
¿Qué valor esperamos que tenga la probabilidad \(P(X\geq 0.7)\)?
Podemos calcularlo con valores lógicos, igual que en el caso discreto
## [1] 0.300253
EJERCICIO 2
Calcula ahora tú:
Seguramente a todos nos suena la campana de Gauss. Esta no es más que la función de densidad de la distribución normal, una de las más utilizadas en estadística. Podemos dibujarla en R de la siguiente manera:
x <- seq(-5, 5, length.out = 10000)
y <- dnorm(x, mean = 0, sd = 1)
plot(x, y, type = "l", col = "blue")Si queremos dibujarla a mano, solo tenemos que tener en cuenta que debemos pintar el máximo de la campana en la media (\(\mu\)), y los puntos de inflexión en la media más menos la desviación típica, \(\mu \pm \sigma\).
Entonces, ¿por qué es tan importante?
Recordemos que la pinta que tenía la variable uniforme continua era la siguiente
Pero, ¿qué pinta tiene la suma de dos uniformes?
Solución
Vamos a volver a verlo numéricamente, primero creamos dos variables, y la suma de las dos
Y luego lo dibujamos
De nuevo es fácil ver en este ejemplo, un triángulo, como el área bajo la función de densidad es uno usando el área de un triángulo.
¿Y de 5?
Solución
Igual que antes, vamos a verlo numéricamente. Creamos 3 variables más, y las sumamos a las dos anteriores
Y luego lo dibujamos
Y vemos que se empieza a parecer a la campana de Gauss. De hecho podemos dibujar la campana encima de la misma y ver cómo de similares son
hist(x_sum_5, freq = FALSE, main = "Suma 5 uniformes")
lines(density(x_sum_5), col = "blue")
x_bar <- mean(x_sum_5)
s <- sd(x_sum_5)
x <- seq(min(x_sum_5), max(x_sum_5), length.out = 10000)
y <- dnorm(x, mean = x_bar, sd = s)
lines(x, y, col = "red")¿Qué es lo que hay detrás?
Solución
Se trata del teorema central del límite. Este nos dice que sumando distribuciones con media (\(\mu\)) y desviación típica (\(\sigma\)), la distribución de la suma se asemeja mucho a una normal.
Podemos ver que sucede con otras distribuciones, como la exponencial. Veamos primero cómo es una exponencial, que tiene la siguiente función de distribución
\[f(x) = \lambda e^{-\lambda x}\]
Y ahora su suma
x_exp2 <- rexp(100000)
x_exp3 <- rexp(100000)
x_exp4 <- rexp(100000)
x_exp5 <- rexp(100000)
x_exp6 <- rexp(100000)
x_exp7 <- rexp(100000)
x_exp8 <- rexp(100000)
x_exp9 <- rexp(100000)
x_exp10 <- rexp(100000)
x_sum <- x_exp1 + x_exp2 + x_exp3 + x_exp4 + x_exp5 + x_exp6 + x_exp7 + x_exp8 + x_exp9 + x_exp10
hist(x_sum, freq = FALSE, main = "Suma 10 exponenciales")
lines(density(x_sum), col = "blue")
x_bar <- mean(x_sum)
s <- sd(x_sum)
x <- seq(min(x_sum), max(x_sum), length.out = 10000)
y <- dnorm(x, mean = x_bar, sd = s)
lines(x, y, col = "red")Dependiendo de la distribución que utilicemos, el número de distribuciones que tenemos que sumar
Podemos ver un ejemplo físico, utilizado para “demostrar” el teorema, conocido como el Galton Board
EJERCICIO 3
Comprobar que ocurre cuando sumamos diferentes variables entre sí: uniforme, exponencial, normal, binomial…
Repetir el proceso anterior e ilustrar gráficamente.
Obtén observaciones al azar de cualquier población de media finita \(\mu\). A medida que el número de observaciones obtenidas aumenta, la media \(\bar{x}\) de los valores observados se acerca más y más a \(\mu\), la media poblacional.
La ley de los grandes números es el principio en el que se basan negocios como los casinos y las compañías de seguros. Las ganancias (o las pérdidas) de un jugador en un juego son inciertas -por este motivo el juego es emocionante- pero las ganancias de muchos jugadores se pueden estimar.
Tengo una población, hago un censo (o tomo una muestra), y espero que la media de la muestra se parezca a mi parámetro. No llego al parámetro, ya que tengo que tomar \(n < N\), el tamaño de la muestra es menor que el de la población. \(n\) depende del dinero/tiempo que tengo.
Ahora, si tomo varias muestras, ¿cuál es el bueno?
La ley de los grandes números garantiza que al medir suficientes objetos, el estadístico se acercará mucho al parámetro desconocido \(\mu\).
Podemos proceder de la siguientes manera:
En la práctica es demasiado caro tomar muchas muestras de poblaciones grandes como, por ejemplo, la de los adultos residentes en la Unión Europea. De todas formas, podemos imitar la obtención de muchas muestras utilizando un programa de ordenador. La utilización de un programa para imitar el comportamiento del azar se llama simulación.
Vamos a ver un ejemplo con el conjunto de datos iris, que son muestras de lirios tomadas por Fisher.
En primer lugar veremos qué pinta tiene el dataframe
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
Y trabajaremos con la variable Sepal.Length
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.300 5.100 5.800 5.843 6.400 7.900
## [1] 0.8280661
Para tomar una muestra, igual que al inicio del tema, podemos usar la función sample:
## [1] 6.3 6.3 6.4 5.7 5.1 6.2 7.7 5.8 6.0
## [1] 5.1 4.9 5.5 6.8 5.7 5.2 5.1 5.6 4.5 4.7 7.7 6.5 6.7 6.7 6.1 4.8 5.2 5.7 7.7
## [20] 5.0 6.7 5.5 5.8 5.6 6.9 6.7 6.8 4.8 5.0 6.7 6.2 5.2 6.2 6.7 5.6 5.5
Ahora, vamos a usar un bucle para imitar la toma de 1000 muestras diferentes
medias <- c()
k <- 1000
n <- 9
for(i in 1:k){
medias <- c(medias,mean(sample(iris$Sepal.Length, n, replace=T)))
}
head(medias)## [1] 6.266667 6.366667 5.700000 5.533333 6.277778 5.644444
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.956 5.644 5.822 5.841 6.044 6.789
## [1] 0.2827933
Y con muestras de tamaño 36
medias36 <- c()
k <- 1000
n <- 36
for(i in 1:k){
medias36 <- c(medias36, mean(sample(iris$Sepal.Length, n, replace=T)))
}
head(medias36)## [1] 5.616667 5.797222 6.038889 6.027778 5.938889 5.875000
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.436 5.753 5.853 5.847 5.944 6.247
## [1] 0.1415552
La media, como vimos, parece que se acerca a la media poblacional, ¿y la desviación típica?
Aquí, al final, estamos aplicando otra vez el Teorema Central del Límite. Eso, en definitiva, nos dice que si la distribución de la población es \(\mathcal{N}(\mu,\sigma)\), la distribución de la media muestral, \(\bar{x}\), tiene una distribución \(\mathcal{N}(\mu,\sigma/\sqrt{n})\), con \(n\) el tamaño de la muestra.
Podemos de hecho ver que realmente \(s\) y \(\sigma/sqrt{n}\) son próximos en los casos que hemos visto antes
## [1] 0.276022
## [1] 0.2827933
## [1] 0.138011
## [1] 0.1415552
EJERCICIO 4
Intenta repetir el proceso para la variable carat del dataframe diamonds y ver cómo, con suficientes muestras, la distribución se aproxima a una normal. Calcula la desviación típica de la muestra, y la que esperarías en función de la desviación de la variable. Prueba con tamaños de muestra 144 o 225.
## # A tibble: 53,940 x 10
## carat cut color clarity depth table price x y z
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43
## 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31
## 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31
## 4 0.290 Premium I VS2 62.4 58 334 4.2 4.23 2.63
## 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75
## 6 0.24 Very Good J VVS2 62.8 57 336 3.94 3.96 2.48
## 7 0.24 Very Good I VVS1 62.3 57 336 3.95 3.98 2.47
## 8 0.26 Very Good H SI1 61.9 55 337 4.07 4.11 2.53
## 9 0.22 Fair E VS2 65.1 61 337 3.87 3.78 2.49
## 10 0.23 Very Good H VS1 59.4 61 338 4 4.05 2.39
## # ... with 53,930 more rows
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2000 0.4000 0.7000 0.7979 1.0400 5.0100
Hemos visto que la distribución normal consta de dos parámetros, la media \(\mu\) y la desviación típica \(\sigma\). Hemos visto también con el Teorema Central del Límite distribuciones con diferentes parámetros en función de las distribuciones de las que partíamos.
Podemos ver cómo afecta un cambio de desviación típica a una variable normal
x <- seq(-10, 10, length.out = 10000)
y1 <- dnorm(x, mean = 0, sd = 0.5)
y2 <- dnorm(x, mean = 0, sd = 1)
y3 <- dnorm(x, mean = 0, sd = 2)
y4 <- dnorm(x, mean = 0, sd = 3)
plot(x, y1, type = "l", col = "blue")
lines(x, y2, col = "green")
lines(x, y3, col = "red")
lines(x, y4, col = "purple")
legend("topright", legend=c("sd = 0.5", "sd = 1", "sd = 2", "sd = 3"), col=c("blue", "green", "red", "purple"), lty = 1, inset = 0.02)O un cambio en la media
x <- seq(-10, 10, length.out = 10000)
y1 <- dnorm(x, mean = -5, sd = 1)
y2 <- dnorm(x, mean = 0, sd = 1)
y3 <- dnorm(x, mean = 2, sd = 1)
y4 <- dnorm(x, mean = 5, sd = 1)
plot(x, y1, type = "l", col = "blue")
lines(x, y2, col = "green")
lines(x, y3, col = "red")
lines(x, y4, col = "purple")
legend("topright", legend=c("mean = -5", "mean = 0", "mean = 2", "mean = 5"), col=c("blue", "green", "red", "purple"), lty = 1, inset = 0.02)Vamos a trabajar con dos normales distintas, una \(\mathcal{N}(0, 1)\), la “estándar”, y otra cualquiera como \(\mathcal{N}(5, 2.5)\).
Para una muestra suficientemnete grande de cada una de ellas, vamos a calcular la probabilidad de que los valores se desvíen de la media por debajo de \(-\sigma\), y por encima de \(\sigma\).
n <- 100000
normal_estandar <- rnorm(n, 0, 1)
normal_5_25 <- rnorm(n, 5, 2.5)
sum(normal_estandar < -1)/n## [1] 0.15752
## [1] 0.15967
¿Cómo esperamos que sean los valores mayores a una desviación típica, \(\sigma\)?
## [1] 0.16087
## [1] 0.15584
Vemos que, dado que la distribución normal es simétrica, los valores son muy similares a los previos.
EJERCICIO 5
Calcular las probabilidades a la izquierda y derecha de \(\pm 2\sigma\) y \(\pm 3\sigma\)
Esto nos da lugar a la llamada regla del 68-95-99.7, que nos habla de los valores comprendidos entre \((-\sigma, \sigma)\), \((-2\sigma, 2\sigma)\) y \((-3\sigma, 3\sigma)\).
Podemos calcular la probabilidad teórica con la función qnorm, a la que le damos una probabilidad y nos da el valor que deja a la izquierda dicha área o probabilidad. Podemos además darle la media y desviación típica de la distribución.
## [1] -0.9944579
## [1] -1
## [1] 2.513855
## [1] 2.5
## [1] -3.696121
## [1] -3.7
O darle la vuelta, la función pnorm, dado un valor, nos da la probabilidad a la izquierda del mismo
## [1] 0.1586553
## [1] 0.1586553
## [1] 0.1586553
Por tanto, parece que las normales son en realidad la misma, con un cambio de escala/forma. Entonces, ¿podemos pasar de una a otra?
Habíamos mencionado antes que la normal “estándar” era \(\mathcal{N}(0, 1)\). Estos son los parámetros más habituales, y la que se usa de referencia, por lo tanto podemos “ir y venir” de esta distribución a una con los parámetros arbitrarios que queramos.
En definitiva, como hemos visto antes, considerando la normal \(\mathcal{N}(-3, 0.7)\), querríamos que el -1 en la estándar se correspondiera con el -3.7 en esta, o el -2 con el -4.4, por ejemplo. ¿Cómo lo hacemos?
Hemos visto, en gráficas anteriores, que la media, \(\mu\), simplemente desplaza en el eje horizontal la curva. Por lo tanto podemos restar la media, teniendo, en los casos anteriores, -0.7 y -1.4 respectivamente. Para llegar al valor en la estándar, tras restar la media, nos quedaría entonces dividir por la desviación típica, \(\sigma\). Podemos ver que ahí sí que llegaríamos a -1 y -2 respectivamente.
La fórmula a aplicar, por tanto, es: dado un \(x\), valor en \(\mathcal{N}(\mu, \sigma)\), podemos obtener el valor correspondiente, \(z\), en \(\mathcal{N}(0, 1)\) a través de la siguiente fórmula
\[z = \frac{x-\mu}{\sigma}\]
Este proceso se llama estandarización.
Por ejemplo, el si quisiéramos ver qué valor se corresponde en una normal estándar a 2 en una \(\mathcal{N}(4, 2.5)\) haríamos
## [1] -0.8
Podríamos hacer el proceso opuesto, “desestandarizar”, pasando a una normal de media \(\mu\) y desviación típica \(\sigma\), despejando en la ecuación anterior
\[x = \sigma z +\mu\]
De forma análoga, para ir de -1 en una \(\mathcal{N}(0, 1)\) al valor correspondiente en \(\mathcal{N}(4, 2.5)\) haríamos
## [1] 1.5
EJERCICIO 6
Encuentra los valores correspondientes en una normal de media 0 y desviación típica 1, y calcula la probabilidad de que x valga menos que dichos valores, a
De forma análoga, encuentra los valores correspondientes en una normal de media 3 y desviación típica 1.5, y calcula la probabilidad de que x valga menos que dichos valores, a
Habíamos visto que al calcular la media de una muestra, la distribución de la media \(\bar{x}\), sigue una distribución \(\mathcal{N}(\mu,\sigma/\sqrt{n})\). Normalmente, \(\mu\) es desconocido, ya que es el parámetro que queremos estimar. Podemos conocer, de nuestra muestra, la media \(\bar{x}\) y la desviación típica \(s\).
Sabiendo esto, y dándole una vuelta de tuerca al proceso de estandarización, que en nuestro caso sería
\[ z = \frac{\bar{x} - \mu}{\sigma/\sqrt{n}}\]
Podemos dar un rango de valores, conocido como intervalo de confianza, en el que con una cierta confianza, normalmente el 95%, esperamos que esté el parámetro real, \(\mu\).
Se trata de despejar, en la fórmula anterior, \(\mu\), y escoger el \(z\) que nos dé la confianza deseada. Sería
\[\mu = \bar{x} \pm z\frac{\sigma}{\sqrt{n}}\]
Habría que elegir el \(z\) que nos deje (1 - confianza)/2 a izquierda y derecha. Si queremos una confianza del 95%, querríamos que a izquierda y derecha dejara una probabilidad de 0.025. Vemos que este valor es
## [1] -1.959964
Aproximadamente \(z = 1.96\) para una confianza del 95%.
Vamos ahora a ver un ejemplo, trabajando con el dataframe flights, y la variable dep_delay
## # A tibble: 6 x 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 -1 1004 1022
## 5 2013 1 1 554 600 -6 812 837
## 6 2013 1 1 554 558 -4 740 728
## # ... with 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
Podemos ver qué ocurre tomando muestras de tamaño 36, 225 y 10.000 con los intervalos de confianza
# Para tamaño 36
n <- 36
muestra1 <- sample(flights$dep_delay, n)
(x_bar1 <- mean(muestra1, na.rm = TRUE))## [1] 10.51429
## [1] 29.75327
## [1] -2.620759
## [1] 23.64933
# Para tamaño 225
n <- 225
muestra2 <- sample(flights$dep_delay, n)
(x_bar2 <- mean(muestra2, na.rm = TRUE))## [1] 8.529148
## [1] 28.76791
## [1] 3.27513
## [1] 13.78317
# Para tamaño 10.000
n <- 10000
muestra3 <- sample(flights$dep_delay, n)
(x_bar3 <- mean(muestra3, na.rm = TRUE))## [1] 12.67135
## [1] 40.58252
## [1] 11.88325
## [1] 13.45945
Vemos que el tamaño de la muestra, que depende del tiempo o dinero que invirtamos en el experimento u observación, afecta al tamaño del intervalo de confianza. De hecho podemos calcular la amplitud del intervalo de confianza antes de tomar una muestra siempre que conozcamos \(\sigma\).
Podemos, además, comprobar que se cumple efectivamente la confianza que hemos escogido. Para ello, vamos a calcular muchos intervalos de confianza, y ver qué porcentaje de veces la media real, \(\mu\), está fuera de ellos.
Para una confianza del 95%
medias_flights <- c()
k <- 1000
n <- 225
for(i in 1:k){
medias_flights <- c(medias_flights,mean(sample(flights$dep_delay, n), na.rm = TRUE))
}
extremo_sup <- medias_flights + qnorm(0.025)*sigma/sqrt(n)
extremo_inf <- medias_flights - qnorm(0.025)*sigma/sqrt(n)
sum(mu < extremo_sup | mu > extremo_inf)/k## [1] 0.058
Y del 90%
medias_flights <- c()
k <- 1000
n <- 225
for(i in 1:k){
medias_flights <- c(medias_flights,mean(sample(flights$dep_delay, n), na.rm = TRUE))
}
extremo_sup <- medias_flights + qnorm(0.05)*sigma/sqrt(n)
extremo_inf <- medias_flights - qnorm(0.05)*sigma/sqrt(n)
sum(mu < extremo_sup | mu > extremo_inf)/k## [1] 0.112
EJERCICIO 7
Usando la variable carat del dataframe diamonds, como en el ejercicio 4, calcula interavalos de confianza para la media muestral tomando muestras de tamaño 9, 144 y 40.000
¿Qué pasa si tampoco conocemos la varianza, \(\sigma\)? En el caso de que no conozcamos la verdadera varianza, hay un resultado que nos dice que
\[ t_{n-1} = \frac{\bar{x} - \mu}{\sigma/\sqrt{n}}\]
Donde \(t_{n-1}\), en lugar de una normal sigue una distribución t de Student con \(n-1\) grados de libertad.
A efectos prácticos, esto simplemente significa que en lugar de calcular \(z\) con qnorm, debemos hacerlo con qt, al que damos también el parámetro de sus grados de libertad.
Como en el ejemplo anterior
# Para tamaño 36
n <- 36
muestra1 <- sample(flights$dep_delay, n)
(x_bar1 <- mean(muestra1, na.rm = TRUE))## [1] 9.5
## [1] 24.5735
## [1] 1.185522
## [1] 17.81448
# Para tamaño 225
n <- 225
muestra2 <- sample(flights$dep_delay, n)
(x_bar2 <- mean(muestra2, na.rm = TRUE))## [1] 10.58371
## [1] 49.4164
## [1] 4.091677
## [1] 17.07574
# Para tamaño 10.000
n <- 10000
muestra3 <- sample(flights$dep_delay, n)
(x_bar3 <- mean(muestra3, na.rm = TRUE))## [1] 12.45378
## [1] 39.228
## [1] 11.68483
## [1] 13.22273
Gráficamente, la t de Student, dependiendo de sus grados de libertad, tiene la siguiente forma
No es difícil intuir que, con infinitos grados de libertad, la t de student coincide con la normal.
Esto se puede interpretar como el hecho de que si yo selecciono toda la población (tamaño de muestra \(n = \inf\)), entonces la desviación típica de la muestra es la misma que la de la media.
EJERCICIO 8
Repite el ejercicio 7, calculando ahora intervalos de confianza considerando la desviación típica como desconocida, y tomando \(s\), la desviación típica de la muestra.